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ABSTRACT 


Accurate  seismic  event  location  is  key  to  monitoring  the  Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT)  and  is 
largely  dependent  on  our  understanding  of  the  crust  and  mantle  velocity  structure.  This  is  particularly  challenging  in 
aseismic  regions,  devoid  of  calibration  data,  which  leads  us  to  rely  on  a  priori  constraints  on  the  velocities.  We 
investigate  our  ability  to  improve  seismic  event  location  in  the  Middle  East,  North  Africa,  and  the  Former  Soviet 
Union  (ME/NA/FSU)  by  using  a  priori  three-dimensional  (3-D)  velocity  models  in  lieu  of  more  commonly  used  one¬ 
dimensional  (1-D)  models.  Event  locations  based  on  1-D  models  are  often  biased,  as  they  do  not  account  for 
significant  travel-time  variations  that  result  from  heterogeneous  crust  and  mantle;  it  follows  that  3-D  velocity  models 
have  the  potential  to  reduce  this  bias.  Here,  we  develop  a  composite  3-D  model  for  the  ME/NA/FSU  regions.  This 
fully  3-D  model  is  an  amalgamation  of  studies  ranging  from  seismic  reflection  to  geophysical  analogy.  Our  a  priori 
model  specifies  geographic  boundaries  and  velocity  structures  based  on  geology,  tectonics,  and  seismicity  and 
information  taken  from  published  literature,  namely  a  global  sediment  thickness  map  of  1°  resolution  (Laske  and 
Masters,  1997),  a  regionalized  crustal  model  based  on  geology  and  tectonics  (Sweeney  and  Walter,  1998; 
Bhattacharyya  et  ah,  2000;  Walter  et  ah,  2000),  and  regionalized  upper  mantle  (RUM)  models  developed  from 
teleseismic  travel  times  (Gudmundsson  and  Sambridge,  1998).  The  components  of  this  model  were  chosen  for  the 
complementary  structures  they  provide.  The  1°  sediment  map  and  regionalized  crustal  model  provide  detailed 
structures  and  boundaries  not  available  in  the  more  coarse  5°  models  used  for  global-scale  studies.  The  RUM  models 
offer  improved  resolution  over  global  tomography,  most  notably  above  depths  of  300  km  where  heterogeneity  is 
greatest;  however,  we  plan  to  test  other  published  upper  mantle  models  of  both  P-  and  5-wave  velocity. 

We  compute  travel  times  through  this  integrated  model  for  comparison  with  other  standard  1-D  models,  as  our  goal  is 
to  evaluate  whether  the  3-D  model  can  better  predict  the  observed  travel  times.  The  arrival  times  are  computed 
through  the  model  using  a  3-D  finite-difference  technique  and  are  then  compared  with  a  declustered  set  of  ISC  P 
arrival  times  (Engdahl  et  ah,  1998).  Our  ME/NA/FSU  model  predicts  the  P  and  Pn  travel  times  very  well,  as 
measured  by  variance  reduction,  for  three  stations  we  tested:  ARU,  KVT,  and  GAR;  these  predicted  times  also 
resemble  some  patterns  seen  in  Pn  tomography  models  of  this  region.  Such  tests  will  allow  us  to  identify  parts  of  the 
model  that  may  need  modification. 

We  also  compute  model-based  correction  surfaces  for  each  station  in  the  ME/NA/FSU  regions  that  can  be  used  as 
additional  constraints  in  our  event  location  algorithm  to  determine  the  improvement  provided  by  using  3-D  models. 
We  use  this  method  to  relocate  a  set  of  ground  truth  events:  the  1991  Racha  aftershock  sequence  which  was 
investigated  by  Myers  and  Schultz  (2000)  using  empirical  kriged  correction  surfaces  and  a  1-D  velocity  model.  They 
find  an  epicenter  mislocation  bias  of  42  km  when  no  corrections  are  applied  and  that  this  mislocation  is  reduced  to 
13  km  when  their  empirically  derived  correction  surfaces  are  included.  We  relocate  this  same  set  of  events  using  our 
model-based  correction  surfaces  and  produce  a  mislocation  bias  of  only  26  km,  a  significant  improvement.  We  are 
currently  implementing  methods  to  quantify  uncertainties  on  the  model-based  corrections  which  will  be  required  to 
compute  representative  error  ellipses  for  the  new  locations.  We  also  plan  to  combine  both  the  model-based  and 
empirical  correction  techniques  to  achieve  the  best  improvement  in  location.  This  test  case  demonstrates  the  power  of 
using  3-D  velocity  models  to  improve  location  capability  for  small,  regionally  recorded  events.  This  example  also 
shows  how  the  model-based  approach  holds  great  potential  for  improving  locations  in  aseismic  regions  where  it  may 
not  be  possible  to  compute  empirical  correction  surfaces. 
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OBJECTIVE 


High  quality  regional  velocity  models  are  necessary  for  accurate  seismic  locations  required  for  monitoring  the 
Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT).  A  common  approach  to  seismic  location  has  been  to  use  a  one¬ 
dimensional  (1-D)  velocity  model  and  apply  source-specific  station  corrections  to  remove  bias  due  to  three- 
dimensional  (3-D)  Earth  structure.  Such  corrections  yield  more  accurate  travel-time  predictions  and  hence  more 
accurate  location  of  the  seismic  event;  however,  they  do  not  provide  source-specific  uncertainties  for  the  new 
locations. 

Empirical  kriged  corrections  and  other  model  optimization  methods  are  well  suited  to  improve  travel-time  prediction 
in  areas  where  ground  truth  (GT)  events  are  available;  however,  vast  portions  of  the  ME/NA/FSU  are  devoid  of  GT 
events  ( e.g .,  North  Africa,  Russia).  In  order  to  improve  travel-time  prediction  at  regional  and  near  teleseismic 
distances  in  regions  without  GT  calibration,  we  have  developed  the  a  priori  3-D  velocity  model  MENAFSU1.0.  We 
seek  to  validate  our  a  priori  3-D  velocity  model  MENAFSU1.0  and  use  this  model  to  locate  a  regionally  recorded 
set  of  GT  events.  The  new  locations  are  compared  to  those  computed  using  a  1-D  velocity  model  and  empirically 
derived  kriged  correction  surfaces,  and  the  improvement  in  location  using  our  model-based  correction  surfaces  is 
demonstrated  in  an  end-to-end  case  study. 


RESEARCH  ACCOMPLISHED 


Middle  East,  North  Africa,  and  Former  Soviet  Union  3-D  velocity  model  (MENAFSU1.0) 

This  regionalized  model  is  a  preliminary  set  of  geophysically  distinct  regions  that  can  be  used  for  estimating  travel- 
times,  surface  wave  dispersion,  and  discrimination  properties  particularly  in  aseismic  regions  where  calibration  data 
is  sparse.  The  model  can  also  provide  a  platform  for  assessing  progress  in  seismic  location,  discrimination, 
detection — the  entire  calibration  process — and  aid  in  determining  the  priority  and  planning  of  calibration 
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Figure  2h  This  volume  of  travel  l  imes  is  Ihc  output  finm  the  T-D  Unite  difference  enBoilnlimi  The  3-D  vekiehy 
structure  causes  live  travel  limes  Id  be  asymmetric  about  thar  Nation  KVT.  which  is  l^ated  al  .zero  lime  t  bright  blue 
spot  in  canter).  Arrival  times  fin  P  and  Pn  phases  far  speeiik  aouree-staakin  paths  are  extracted  from,  this  volume  by 
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experiments.  Our  goal  is  to  build  a  reference  model  for  the  region  that  can  be  used  in  making  model-based 
correction  surfaces.  This  model  has  been  used  to  derive  predictions  of  surface  wave  dispersion  curves  and  regional 
travel  times  for  comparison  with  data.  In  addition,  empirical  observations  are  being  used  to  further  refine  the 
background  geophysical  model.  The  MENAFSU1.0  model  specifies  geographic  boundaries  and  velocity  structures 
based  on  geology,  tectonics,  and  seismicity.  This  regionalization  serves  as  a  starting  point,  and  we  expect  to  refine 
and  improve  upon  it  based  on  such  tests  as  predicting  P  and  Pn  travel  times  (presented  here)  as  well  as  surface  wave 
velocities.  As  the  model  improves  and  demonstrates  some  predictive  power,  it  may  evolve  into  a  base  model  for 
tomographic  inversions  as  well  as  a  reference  model  for  other  CTBT-related  research  efforts. 

The  entire  area  is  divided  into  3 1  regions  selected  primarily  on  the  basis  of  seismicity  patterns,  tectonic  history, 
crustal  thickness,  and  crustal  velocity  (Bhattacharyya  et  ah,  2000;  Walter  et  ah,  2000).  Each  region  is  composed  of 
7  homogeneous  layers  for  water,  unconsolidated  and  consolidated  sediments,  upper  crust,  middle  crust,  lower  crust, 
and  upper  mantle  half-space,  following  the  specifications  of  Mooney  et  ah  (1998)  in  their  5°  resolution  model 
CRUST5.1.  We  have  stripped  off  and  replaced  the  sediments  with  a  high-resolution  1°  sediment  model  (Laske  and 
Masters,  1997)  and  have  discretized  all  31  regions  into  1°  latitude  and  longitude  blocks  while  the  depth  layers 
remain  variable  in  thickness.  A  representative  cross-section  of  the  velocity  model  is  shown  in  Figure  1.  We  have 
added  a  regionalized  upper  mantle  model  below  the  Pn,  Sn  velocity  layers.  Of  the  many  upper  mantle  models 
available,  we  have  chosen  initially  to  use  the  RUM  (Regionalized  Upper  Mantle)  model  of  Gudmundsson  and 
Sambridge  (1998).  This  model  is  attractive  for  two  reasons;  first  it  follows  a  similar  regionalization  procedure  as 
was  used  for  the  crust  except  applied  to  the  mantle.  Second,  it  is  based  on  deviations  from  AK135  (Kennett  et  ah, 
1995)  which  is  our  basis  for  calculating  travel  times. 

The  MENAFSU1.0  velocity  model  is  an  amalgamation  of  geophysical  and  geological  studies,  constituting  a  self- 
consistent  3-dimensional  Earth  model  for  the  crust  and  upper  mantle.  Because  the  MENAFSU  model  is  3-D, 
region-specific  velocity  structure,  it  can  be  characterized  more  accurately  than  1-  or  2-D  models.  While  this  model 
comprises  P- wave,  5-wave,  density,  and  attenuation,  we  use  only  the  P- wave  structures  at  this  time. 

Finite  Difference  Travel-Time  Calculations 

We  use  an  algorithm  originally  developed  by  Vidale  (1988)  and  further  refined  by  Hole  and  Zelt  (1995)  which  uses 
a  finite  difference  (FD)  approximation  to  the  Eikonal  equation  to  compute  first  arrival  travel  times  through  regularly 
gridded  velocity  structures.  This  technique  approximates  Huygen's  principle  by  propagating  wavefronts  radially 
outward  from  a  point  source  using  each  grid  (or  time)  node  as  a  secondary  source  for  each  successive  grid  node. 

This  procedure  is  more  efficient  and  accurate  than  ray  tracing  as  it  is  able  to  account  for  sharp  velocity  gradients 
which  can  produce  different  propagation  modes  (e.g.,  refracted  and  diffracted  body  waves,  head  waves)  in  addition 
to  direct  phases.  Furthermore,  because  the  travel  times  are  computed  for  every  grid  point  in  the  volume,  the  code  is 
much  faster  than  tracing  rays  from  a  single  event  to  a  large  number  of  receivers.  In  fact,  we  use  the  principle  of 
reciprocity  by  placing  the  source  point  at  the  seismic  station  location  and  interpolate  the  travel  times  computed  at 
each  grid  node  to  match  the  exact  earthquake  location  (Flanagan  et  ah,  1999). 

We  modify  the  original  code  in  two  ways.  First,  we  adapt  it  to  read  in  3-D  velocity  models  instead  of  1-D  such  that 
it  can  compute  times  through  our  MENAFSU  1.0  (or  any  custom  3-D)  model.  Second,  we  apply  a  Cartesian  to 
spherical  coordinate  transformation  to  the  source  and  receiver  locations  that  are  input  to  the  code  (K.  Koper, 
personal  communication,  1999).  These  modifications  are  necessary  as  we  need  to  compute  travel  times  out  to 
regional  and  near-teleseismic  distances  (-13°  to  30°)  while  the  original  FD  algorithm  operates  in  Cartesian 
coordinates,  which  are  accurate  only  for  distances  of  up  to  300  km  in  the  Earth.  To  account  for  the  curvature  of  the 
Earth  one  option  is  to  apply  an  Earth  flattening  transformation;  however,  this  is  only  strictly  valid  for  1-D  velocity 
models  and  not  for  our  3-D  model.  We  develop  the  spherical  coordinate  transformation  scheme  in  which  we 
parameterize  a  spherical  surface  inside  the  Cartesian  volume  of  grid  nodes.  Here,  the  source  and  receiver  positions 
and  velocity  model  are  known  in  spherical  coordinates,  and  they  are  transformed  into  the  Cartesian  system  as  input 
to  the  FD  code.  The  finite  difference  operators  remain  unchanged,  and  the  output  is  converted  to  spherical 
coordinates  as  well. 

The  FD  algorithm  requires  our  velocity  model  to  be  discretized  on  to  an  equally  spaced  grid;  this  requires  fine 
sampling,  by  linear  interpolation,  of  our  1°  by  1°  velocity  model.  This  procedure  results  in  the  loss  of  some 
resolution  with  depth  (e.g.,  some  sediment  layers  may  be  thinner  than  the  grid  size)  and  an  oversampling  of  lateral 
structure  as  the  velocities  in  the  shallower  part  of  the  model  vary  more  rapidly  with  depth  than  they  do  laterally. 


Our  fully  3-D  code  is  run  in  a  volume  of  dimensions  of  roughly  30°  by  30°  laterally  and  800  to  1000  km  deep  with  a 
grid  spacing  of  3  km;  the  exact  latitude -longitude  bounds  of  the  volume  vary  with  each  station.  The  grid  spacing  is 
determined  empirically  as  a  trade-off  between  the  accuracy  of  the  travel-time  prediction  and  computer  memory 
limitations  as  both  the  velocity  grid  (or  slowness)  and  the  travel-time  grid  must  be  held  in  memory  simultaneously. 
Our  current  resources  limit  us  to  grids  with  about  350  million  nodes  which  requires  approximately  4  Gbytes  of 
computer  memory.  After  a  number  of  tests  running  the  code  with  incrementally  smaller  grid  sizes  and  comparing 
the  resulting  travel  times  with  those  computed  using  a  simple  1-D  ray  tracer,  we  find  that  a  grid  spacing  of  3  km 
provides  a  reasonable  accuracy  as  shown  in  Figure  2a.  Accuracy  and  errors  of  the  FD  code  arise  from  grid  boundary 
effects  and  from  the  discrete  approximation  of  the  differencing  algorithm.  These  accuracies  are  being  evaluated 
through  convergence  tests  and  through  careful  comparisons  with  known  solutions;  currently  the  RMS  errors  appear 
to  be  less  than  0.5  seconds  (Figure  2a).  Absolute  errors  and  accuracy  are  being  derived  in  seismic  areas  where  we 
can  statistically  assess  the  ability  of  this  a  priori  model  to  predict  data. 

The  outputs  from  the  FD  code  are  two  3-D  volumes:  one  of  the  P-wave  velocity  at  each  grid  node  and  one  of  the  P 
arrival  time  at  each  node;  an  example  of  the  travel-time  volume  is  shown  in  Figure  2b  as  computed  for  station  KVT. 
We  perform  two  runs  for  each  station  of  interest,  once  using  the  MENAFSU1.0  and  once  using  the  global  reference 
model  AK135\  we  then  compare  times  from  each  of  these  volumes  for  the  same  source-station  paths. 

Validation  of  Travel-Time  Improvement 

To  test  the  predictive  power  of  the  MENAFSU1.0  model,  we  compute  the  median  residuals  between  the  observed  P 
and  Pa  arrivals  at  each  station  and  the  arrivals  predicted  by  both  our  MENAFSU1.0  model  and  the  AK135.  We  use 
P  and  Pn  travel-time  picks  from  a  relocated  subset  of  earthquakes  in  the  ISC  Bulletin  (Engdahl  et  al.,  1998)  which 
have  been  further  processed  with  a  declustering  algorithm.  This  entails  filtering  spatial  data  with  the  goal  of 
consolidating  heavily  clustered  points  into  small  collections  of  representative  data  points  with  associated 
uncertainties.  This  technique  is  intended  to  reduce  the  extremely  large  set  of  arrival  times  into  a  smaller  more 
manageable  data  set  with  minimal  loss  of  information.  As  part  of  the  process,  outliers  within  data  clusters  are 
removed  and  the  uncertainty  estimates  for  the  improved  data  set  are  provided. 

Travel-time  residuals  are  computed  by  interpolating  between  grid  nodes  to  calculate  the  predicted  arrival  time  for  an 
exact  source-station  path.  This  is  done  for  both  the  MENAFSU1.0  travel-time  volume  and  the  AK135  travel-time 
volume,  and  each  predicted  time  is  subtracted  from  the  observed  arrival  time.  These  residuals  (observed  -  predicted 
time)  for  station  ARU  and  KVT  are  shown  in  histogram  form  in  Figure  3.  The  3-D  MENAFSU1.0  model  predicts 
the  observations  very  well,  showing  a  significant  variance  reduction,  while  for  others  it  does  not  improve  the  fit 
given  by  the  AK135.  This  procedure  will  help  us  identify  regions  where  the  model  needs  to  be  improved 
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Model-Based  Correction  Surfaces 

To  compute  model-based  correction  surfaces  we  subtract  the  AK135  predicted  time  from  the  MENAFSU  predicted 
time  along  a  regular  grid  in  latitude,  longitude  and  at  a  depth  of  10  km.  Examples  of  these  surfaces  are  shown  in 
Figure  4  for  the  six  stations  we  will  use  in  our  relocation  test  described  below.  Each  correction  surface  provides 
station-specific  travel-time  corrections  for  regional  to  near-teleseismic  distances  providing  improved  accuracy  and 
precision  for  travel-time  predictions.  We  find  travel  time  differences  of  up  to  6  sec  relative  to  AK1 35,  most  in  areas 
of  very  thick  crust  or  sediment  (blue  indicates  fast  regions  and  red  indicates  slow).  Note  the  patterns  in  these 
correction  surfaces  correlate  with  the  structural  features  in  the  MENAFSU1.0  model;  fast  predictions  are  seen  to  the 
north  on  the  Russian  platform  while  slow  anomalies  are  seen  at  the  southern  Caspian  as  well  as  eastern  Turkey  and 
the  southern  Caucasus. 

Next  we  determine  the  improvement  provided  by  the  3-D  models  by  relocating  a  set  of  ground  truth  events.  We  use 
the  model-based  correction  surfaces  computed  for  six  stations  (KAS,  KVT,  GAR,  KHO,  ARU,  and  SVE) 


comprising  the  test  network  which  recorded  the  1991  Racha  aftershock  sequence;  these  surfaces  are  plotted  in 
Figure  4.  The  fast  and  slow  residual  patterns  seen  in  these  correction  surfaces  show  strong  similarities  to  the 
patterns  seen  in  the  empirical  kriged  correction  surfaces  based  on  1-D  velocity  models  Myers  and  Schultz  (2000). 

Improving  Location:  The  1991  Racha  Aftershock  Sequence 

The  improvement  in  seismic  location  that  is  gained  by  using  our  a  priori  3-D  velocity  model  is  the  ultimate  test  of 
this  method's  utility.  We  use  a  set  of  GT2  locations  determined  using  regionally  recorded  aftershocks  in  the  region 
of  Racha,  Georgia.  These  GT  locations  are  then  compared  with  our  locations  computed  using  our  model-based 
correction  surfaces,  and  we  assess  the  degree  to  which  our  3-D  velocity  model  accounts  for  some  of  the  differences 
that  are  indicated  by  the  path  corrections  and  spatial  differences  in  arrival  times  to  each  of  the  six  stations 
(Flanagan,et  ah,  2000). 

In  a  previous  study,  Myers  and  Schultz  (2000)  demonstrated  location  improvement  using  Modified  Bayesian 
Kriging  (MBK)  to  compute  empirical  correction  surfaces  for  the  sparse  6-station  test  network.  The  MBK  correction 
surfaces  for  the  test-network  stations  are  based  on  high-quality  teleseismically  constrained  hypocenters  throughout 
the  Middle  East  (Engdahl  et  ah,  1998).  The  1991  Racha  events  are  then  relocated  with  and  without  the  aid  of  MBK 
correction  surfaces,  and  the  resulting  epicenters  are  compared  to  the  benchmark  GT2  locations  determined  from  a 
dense  local  deployment  of  seismic  sensors.  When  no  travel-time  correction  is  applied,  the  mean  horizontal  distance 
between  the  local  and  test  network  locations  is  42  km,  and  there  is  a  distinct  bias  in  sparse-network  locations 
towards  the  north-northwest.  The  mean  difference  between  local  and  sparse  network  locations  is  reduced  to  13  km 
when  the  empirical  corrections  are  applied  and  the  bias  in  location  is  significantly  reduced  (Figure  5,  center). 


Test  Case:  Rdocalron  of  1991  Racha  earthquake  sequence 
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Next  we  relocate  the  Racha  events  with  and  without  our  new  model-based  correction  surfaces  (shown  in  Figure  4)  as 
constraints,  and  the  resulting  epicenters  are  compared  to  benchmark  GT2  locations.  The  mean  difference  between 
the  benchmark  and  test  network  locations  is  reduced  from  42  km  to  26.6  km  when  the  model-based  corrections  are 
applied,  and  the  bias  in  location  is  significantly  reduced  (Figure  5,  right).  This  reduction  is  almost  as  large  as  that 
achieved  when  the  empirical  kriged  correction  surfaces  are  applied.  We  are  currently  implementing  methods  to 
derive  uncertainty  estimates  for  the  model-based  correction  surfaces  which  will  be  required  to  compute 
representative  error  ellipses  for  the  new  locations.  Since  our  a  priori  MENAFSU1.0  model  does  not  have 


uncertainties  associated  with  it,  we  must  use  the  misfit  between  the  observed  travel  times  and  those  predicted  from 
our  model  as  a  measure  of  uncertainty  (see  Figure  3). 


This  test  case  demonstrates  the  power  of  3-D  model-based  corrections  to  improve  location  of  small,  regionally 
recorded  events.  This  model-based  approach  shows  enormous  potential  for  improving  location  in  aseismic  regions, 
where  GT  events  are  not  available.  However,  this  relocation  of  the  Racha  sequence  should  be  viewed  with  caution 
as  we  continue  to  test  and  evaluate  the  predictive  power  of  the  MENAFSU1.0  velocity  model  in  other  regions. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Empirically  derived  correction  surfaces  are  applicable  where  GT  events  exist,  but  extrapolation  of  travel-time 
corrections  from  GT  events  is  quickly  damped  to  the  background  statistics  of  the  AK135  model  in  regions  where 
there  are  no  data  to  constrain  the  corrections.  Therefore,  model-based  corrections,  which  are  shown  here,  are  better 
suited  to  improve  travel-time  prediction  and  thus  location  capability  in  areas  devoid  of  GT  events. 

The  test  case  relocation  of  the  1991  Racha  sequence  should  be  viewed  as  a  one -point  validation  of  our  3-D  model. 
Additional  validation  of  the  predictive  power  of  the  MENAFSU 1 .0  velocities  is  being  actively  pursued  in  other 
high-quality  GT  regions.  Further  testing  of  the  FD  code  for  grid  size  vs.  accuracy  of  the  predicted  travel  time  is 
underway  as  we  investigate  variable  grid  spacing  techniques  to  allow  us  to  compute  times  for  smaller  grid  spacing 
and  geographic  areas  larger  than  30°.  We  are  also  compiling  a  large  data  set  of  shear  wave  arrival  times  to  test  and 
validate  the  .S'-wave  velocity  structure  of  MENAFSU1.0  as  our  FD  code  can  easily  compute  first  arrival  times  for 
shear  waves. 

This  work  was  performed  under  the  auspices  of  the  U.S.  Department  of  Energy  by  University  of  California 
Lawrence  Livermore  National  Laboratory  under  contract  No.  W-7405-Eng-48  for  the  Office  of  Research  and 
Development,  NN-20,  within  the  Office  of  Nonproliferation  and  National  Security,  NN-1.  LLNL  contribution  UCRL- 
JC-1 38984. 
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